ARPES Signatures of Few-Layer Twistronic Graphenes

Diverse emergent correlated electron phenomena have been observed in twisted-graphene layers. Many electronic structure predictions have been reported exploring this new field, but with few momentum-resolved electronic structure measurements to test them. We use angle-resolved photoemission spectroscopy to study the twist-dependent (1° < θ < 8°) band structure of twisted-bilayer, monolayer-on-bilayer, and double-bilayer graphene (tDBG). Direct comparison is made between experiment and theory, using a hybrid k·p model for interlayer coupling. Quantitative agreement is found across twist angles, stacking geometries, and back-gate voltages, validating the models and revealing field-induced gaps in twisted graphenes. However, for tDBG at θ = 1.5 ± 0.2°, close to the magic angle θ = 1.3°, a flat band is found near the Fermi level with measured bandwidth Ew = 31 ± 5 meV. An analysis of the gap between the flat band and the next valence band shows deviations between experiment (Δh = 46 ± 5 meV) and theory (Δh = 5 meV), indicative of lattice relaxation in this regime.


Sample fabrication
Samples were fabricated in an argon atmosphere using a modified PMMA-based tear-and-stack technique [1], controlled by a remote micromanipulation rig. Pre-prepared hBN on graphite flakes (on SiO 2 substrates) were used as an adhesive layer to tear the graphene flakes supported on a PMMA membrane, allowing accurate control over the twist angle. Samples were annealed in UHV at 300 • C for several hours prior to measurement. For further details see SI section 2.
Angle resolved photoemission spectroscopy ARPES experiments were performed at the nanoARPES branch of the I05 beamline of Diamond Light Source. A choice of two focusing optics are available to perform spatially resolved ARPES: a Fresnel zone plate for submicrometre spatial resolution, and a capillary mirror for improved flux and energy resolution (∼4 µm spatial resolution). All photoemission spectra in the main text were measured using the capillary mirror, apart from the tMBG data in the second panel of Fig. 2e which was measured with the zone plate, with a 90 eV photon energy and linearly polarised light, at a measured sample temperature of < 85 K. Experimental constant energy cuts were averaged over ±5 meV of the stated energy. For further details see SI section 3.
Tight-binding modelling of electronic structure A hybrid k · p theory-tight-binding Hamiltonian was used for the twisted structures, as previously reported [2][3][4]. The SWMcC parameters used for this description are shown in Table 1 and are taken from [5]. Further details are given in SI section 4.
ARPES simulations. The general form of the ARPES intensity in the central mBZ is written as: where H is the twisted graphene Hamiltonian [6] (see SI section 4).The initial state, ψ band , is the wavefunction of the graphene electron, written as a linear combination of Bloch functions at a given point in momentum space, comprised of layer and sublattice components that are coupled using SWMcC parameters as well as mixed by the moiré [7]. The components are solved by diagonalising the respective system Hamiltonian and solving for its wavefunctions. From this, the weights (c l,λ ) are computed for each lattice site, where l is the layer number and λ is the sub-lattice index. The final state, ψ vac , is assumed to be a plane wave in the vacuum. The interaction term ∇ k H adds a small phase shift to the ARPES spectra. Attenuation and interference after photoemission are accounted for by the term: F = Ae ikz·d , where A = 0.4 per graphene layer was determined by comparison to experiment. d = 3.35Å is the distance between adjacent layers and k z is the out-of-plane component of the final state momentum, calculated by determining k z from the conservation of energy (E p + W = ϵ q −hω) and the kinetic energy of an emitted photoelectron (E p =h 2 (k 2 z + k 2 ∥ )/2m e ). The Lorentzian factor L(E p + W − ϵ q −hω) broadens the spectra to match the experimental broadening of 60 meV. Further details are given in SI section 5.  Figure S1: Sample summary. a Low magnification optical microscope image of a sample with a single, grounding contact (ungated sample Sample fabrication was performed using a remotely controlled micromanipulation rig housed inside an argon atmosphere. Twisted multilayers of graphene (including tBG, tMBG and tDBG) were transferred onto hBN/Si(SiO 2 -290nm) and hBN/graphite/Si(SiO 2 -290nm) stacks for ungated and gated samples respectively. The hBN/graphite heterostructures were prepared using a standard PMMA-based dry transfer technique [1], including mechanical exfoliation of crystals on silicon coated with a sacrificial poly-vinyl alcohol (PVA) layer as well as a poly-methyl methacrylate (PMMA) carrier layer. For the twisted graphene multilayers, we employed a modified tear-and-stack technique [8] to the PMMA-based dry transfer method to manipulate the twist angle between the graphene layers. The graphene was placed over the edge of the hBN (with the other half touching the SiO 2 /Si and avoiding the graphite for the gated samples) to tear the graphene layer in half. The sample was then rotated by a target twist angle, θ, before stacking the 2nd half of the graphene layer on top of the first. A single flake of monolayer or bilayer graphene was used for tBG and tDBG twisted samples, respectively. For tMBG twisted samples, flakes consisting of a monolayer attached to bilayer regions were used. The top twisted graphene layer (and bottom graphite flake for gated samples) were contacted using Ti (3 nm)/Au (40 nm) electrodes deposited through a TEM grid shadow mask to minimise contamination. Optical images of ungated and gated samples are shown in Fig. S1. Gated samples were mounted into chip carriers (Fig. S1c) using a room temperature curing, UHV compatible, non-conductive epoxy purchased from Atom Adhesives (AA-bond 2116). Electrical connections between the chip carrier and sample electrodes were made using a wire bonder. Samples were transferred to the beamline in air. Prior to measurement, samples were annealed in UHV at 300 • C for several hours to remove surface adsorbates. Samples without gate electrodes were annealed for ∼3 hours, while samples mounted in the chip carriers were annealed for at least 6 hours.

µARPES
ARPES experiments were performed at the nanoARPES branch of the I05 beamline of Diamond Light Source. Here, synchrotron light is focused to a small spot size using specialised optics to allow for high spatial resolution ARPES. A choice of two focusing optics are available: a Fresnel zone plate for submicrometre spatial resolution, and a capillary mirror for improved flux and energy resolution (∼4 µm spatial resolution). All photoemission spectra in the main text were measured using the capillary mirror (apart from the bilayer on monolayer graphene data in the second panel of Fig. 2e with a 90 eV photon energy and linear horizontal polarised light. Data was collected using a Scienta Omicron DA30 hemispherical analyser, with a chamber pressure of 1 x 10 −10 mbar and a measured sample temperature <85 K. The sample was located and mapped in the ARPES chamber using scanning photoemission microscopy (SPEM). In SPEM, the sample was raster scanned under the focused beam whilst collecting an I(E, k) spectrum at each position, where the direction of k is determined by the orientation of the analyser relative to the sample, giving a 4D dataset of I(E, k, x, y). Each pixel of the SPEM image shown in Fig. S1f is an integrated intensity of the I(E, k) spectrum from that point on the sample. By comparing this with optical images collected during fabrication, different regions on the sample can be easily identified. kspace mapping was performed by rotating the chamber, including the analyser, around the fixed sample and optic, to measure from different polar emission angles. Electrostatic gating measurements were performed using a Keithley 2634B Source Meter, applying a voltage to the graphite back gate layer, while a separate ground contact connects to the twisted graphene layers adjacent to the measurement position (e.g. Fig. S1e). Achieving an effective gate in-situ requires good electrical isolation between the two contacts, requiring precise alignment of the contact pads relative to the flakes as well as careful handling. It was found that not all the samples fabricated with gates could be electrically biased in-situ.
Samples of twisted graphene fabricated by exfoliation and mechanical stacking are known to show heterogeneity at the micrometre scale, with domains of different angular orientation separated by wrinkles, bubbles, and, sometimes, regions displaying commensurate to incommensurate transitions [9].
To assess the uniformity of our samples, we conducted AFM measurements and SPEM mapping prior to acquiring detailed spectra from specific positions. This is summarised in Fig. S2 for a tDBG sample. From the AFM images, uniform areas of interest a few micrometres across can be identified that are free from wrinkles and bubbles (Fig. S2b). On some samples, AFM brushing ('nano-squeegee' [10]) was used in specific locations to remove contamination between the layers and further improve the uniformity, as indicated by the arrows in Fig. S2b. SPEM mapping was used to locate these regions for ARPES analysis.
By integrating over a selected area of energy and momentum in the I(E, k) spectra at each point on the sample, Fig. S2h(i) red dashed boxes, different flakes can be highlighted in the SPEM images (Fig. S2d,e). Comparison with optical and AFM images allows determination of the areas of interest. Averaging the I(E, k, x, y) dataset within these positions indeed shows them to be uniform -sharp bands and clear hybridisation gaps are maintained (Fig. S2f,g). To further demonstrate this, a series of ARPES spectra from different positions are shown in Fig. S2h. When in a heterogeneous area, the spectra appear broadened, and contributions from different twist domains gives a superposition of bands with no clear interactions. By applying this simple spatial analysis, we could identify uniform areas of different twist angle within our samples and identify regions of different twist angles on a given sample. The data presented in the main text come from 17 different measurements (I(E, k x , k y ) spectra from different twist-angles and/or stacking geometries) which were acquired from 9 different samples (4 tBG samples, 4 tMBG samples, 1 tDBG sample).

Max
Min Figure

Continuum model for twisted few-layer graphene
For bilayer graphene (BLG), the standard Hamiltonian as defined by McCann et al. [11] was used. For each twisted structure, a hybrid k·p-tight binding Hamiltonian was used for the twisted structures, similar to those previously reported [2][3][4]. Such Hamiltonians consider the dispersion around a given valley in reciprocal space and information can be obtained about the system by solving for the eigenfunctions and eigenvalues of the system. The Hamiltonians for twisted bilayer graphene (tBG), twisted monolayerbilayer graphene (tMBG) and twisted double-bilayer graphene (tDBG), respectively are given below: where p x , p y are the in-plane components of the momentum shifted to be centred around the K + valley (p =hk −hκ + , where κ + = (4π/3a, 0) ). In Eq. (S5), a second-order expansion in momentum is used to account for trigonal warping of the band structure at higher energies. The interlayer coupling matrices T ij across the twisted interface are defined as: , with j = 0, 1, 2 and K = 4π/3a. Table 2 gives the values for the Slonczewski-Weiss-McClure (SWMcC) parameters used for all of the calculations here.
where v i = √ 3aγ i /2h. v accounts for the intralayer coupling between the A and B sites in the graphene lattices. γ 1 corresponds to the interlayer coupling between the dimer sites (A b /B t ) of the aligned bilayers and is also used to compute the interlayer coupling at the twisted interfaces. v 3 , v 4 account for the coupling between the non-dimer sites and the coupling between a dimer site and a nondimer site, respectively. Lastly, ∆ ′ accounts for the difference in the on-site electron energies of the dimer sites. v 4 and ∆ ′ also have the effect of adding electron-hole asymmetry into the optical response of the system.

Simulating ARPES intensity
The ARPES intensity is proportional to the square of the modulus of the transition amplitude between the initial and final states of the system under the perturbation caused by incoming photons [12,13]. We define the initial state of the system as a single Bloch electron at the surface. The general form of a Bloch wave is: where j is the band index of the system and C j is the amplitude of the wavefunction of the system, containing a sublattice phase. The wavefunctions, φ k,j (r), are calculated by diagonalising the Hamiltonians shown in Eqs. (S2),(S3) and (S4). G m are the moiré reciprocal lattice vectors. The final state of the system, the emitted photoelectron with momentum p e , is treated as a plane wave with approximate form: The perturbation due to the incoming radiation is treated as a first-order approximation of the standard perturbative Hamiltonian: where A is the electromagnetic vector potential of the incoming radiation, H is the Hamiltonian of the graphene system, and v is the electron velocity operator which introduces the interaction Hamiltonian (H int = ∇ k H) for the photoemission process. Accounting for multi-layer effects, the general form of the ARPES intensity in the central mini-Brillouin zone can then be written as: wherehω is the energy of the photons used in the experiment, W = 4.6 eV is the workfunction of graphene [14], and ϵ M is the energy of the measured electron in the crystal lattice. The delta function comes from treating a single electron in a many-body system using a time-ordered one-electron Green's function [15]. To give a qualitative match to the experimental spectra, accounting for instrument resolution and lifetime broadening etc., a Lorentzian with 60 meV broadening is included (L(ϵ e +W −ϵ M −hω)).
A is a 2x1 vector whose value can be changed to alter the polarisation of the light in the system. Eq. (S10) thus gives an analytical approximation to the ARPES intensity.
Accounting for attenuation and interference across multi-layers We do not attempt a onestep photoemission model [15], but instead account for attenuation and interference of the photo-emitted electrons from different layers using a scaling factor, F , assuming a plane wave final state. This scaling factor accounts for attenuation (due to the short mean free path of photo-electrons) and interference (due to layer-dependent phase differences of the photo-electrons): F = Ae ikz·d , where k z is the outof-plane momentum of the photo-emitted electrons and d = 3.35Å is the distance between adjacent graphene layers. The attenuation factor here was set to be A = 0.4 per layer, found by comparison to the experimental data. Starting with the kinetic energy of the emitted photo-electron using conservation of energy, the out-of-plane momentum is This defines all terms in the factor F l−1 = (Ae ikz·d ) l−1 seen at the end of Eq. (S10), where l is the layer number counting from the top surface. As a consequence of Eq. (S12), the phase difference between photoelectrons emitted from different layers depends on the energy of the incident photon used. This phase factor is plotted as a function of the photon energy in Fig. S3, for photoelectrons emitted from neighbouring graphene layers, across the experimentally relevant photon energies. Note that the phase factor changes significantly with photon energy, indicating that the relative intensity of different features in the ARPES spectra should be strongly dependent on the photon energy. The validity of this approach was tested by comparing experimental and predicted spectra for bilayer graphene and tMBG at varying photon energy, as shown in Fig. S4. While many of the key features are consistent across the different energies, the relative intensities of features changes with photon energy. The results show that although this simple model is not a full description of the photoemission process, it gives a good approximation. We also note that the elongated flat region at the beginning of Fig. S3 corresponds to the fact that it takes ϵ = W to photoexcite electrons from the Γ point but additional energy to photoexcite electrons with nonzero in-plane momentum. We note that the comparison between simulated and experimental spectra in Fig. S4 can also be used to confirm the sign of the tight binding parameter γ 1 . There is no clear consensus on the sign of γ 1 in the literature: in some reports it is assumed to be positive, and in others, negative [12]. In other reports, γ 1 is taken to be negative with a phase factor of e iπ per layer replacing the factor exp(ik z · d) in F [16]. A comparison of simulated and experimental spectra demonstrates that the sign of γ 1 is positive with no additional phase needed. We determined the twist angles directly from the ARPES spectra. Replicas bands at E F are positioned at the mBZ corners, which form a honeycomb pattern with the primary bands (Fig. S5a,b). These periodic features were used to determine the mBZs, which provide the κ 1 and κ 2 positions. The twist angle can then be calculated from (S13) The error on the twist angle was estimated from the uncertainty in positioning the mBZs due to the finite width of primary and replica bands at E F . For the case of 1.5 • tDBG, with the flat band at E F replicas are not readily resolved at the mBZ corners. Instead, we used the periodicity of the intensity at γ to position the mBZ, both in the constant energy cuts (Fig. S5e) and the dispersion of the flat band (Fig. S5f). Again, the error bar on this comes from the uncertainty in positioning the mBZs.

Associating spectral features with electronic bands
For small-twist angle samples at higher binding energies, the spectra become more complex, as shown in Figs. S6a,b, and it becomes difficult to distinguish which layer and band each spectral feature is associated with. Comparison between the simulated and experimental spectra, alongside the band structure predictions, enables the ARPES spectral features to be assigned to distinct valence bands. For example, in Fig. S6, constant energy cuts of the ARPES spectra are shown alongside the band structure in the k x − k y plane averaged over the same energy window, E ± 30 meV, shown in the extended zone scheme such that the pattern is repeated across successive mBZs. In Figs. S6e,f the bands that are most apparent in the ARPES spectra are highlighted. In the band calculation the colour scheme is as follows: blue, green, red, cyan = 1st, 2nd, 3rd, 4th valence band.

Analysis of replica band intensity
As shown in Section 7, in the extended zone scheme, bands are equivalent in each mBZ. However, the photoemission intensity of the replica bands decreases in successive mBZs due to the reduced probability of scattering to higher wave vectors. Fig. S7a shows the photoemission intensity near the Fermi level for a 3.0 • tBG. The primary Dirac points are labelled following the convention from the main text of κ 1 for the top layer and κ 2 for the bottom layer, where these now refer to an intensity. The photoemission intensity from the bottom layer is reduced relative to the top layer due to attenuation, i.e. κ 2 = A × κ 1 , where A is the previously mentioned attenuation set to match the experiment as A = 0.4. Intensity due to replica bands can be seen at the corners of the mBZs. These are labelled by the order of their respective intensity, i.e. R 1 is the most intense replica and R 14 is the least intense. This is shown more clearly in Fig. S7b, where the total intensity is taken from the circular regions marked in Fig. S7a, and normalised with respect to the intensity at κ 1 . This is compared to the experimental results where there is relatively good agreement and they roughly follow the same hierarchy.
The intensity ordering of the replicas is nontrivial. For example, we can associate replicas R 1 and R 2 with states from the bottom layer, and replicas R 3 and R 4 with states from the top layer. We know this because, in the case of tMBG, R 1 and R 2 show a bilayer graphene-like dispersion, while R 3 and R 4 show a monolayer graphene-like dispersion. Despite this, R 1 and R 2 show a greater intensity than R 3 and R 4 , even though states belonging to the bottom layer would be expected to be attenuated. In addition to this, replicas which we would expect to be equivalent with respect to their scattering distance from their respective primary show differing intensity. This is seen most clearly for the replicas R 1 , R 2 , R 7 , R 9 , R 12 and R 14 . Although they are all a single moiré reciprocal lattice vector away from κ 2 , R 7 , R 9 , R 12 and R 14 are all at least an order of magnitude weaker than R 1 and R 2 . The good agreement between experiment and theory confirms the validity of the model used for simulating the photoemission intensity.   3.0° (th) 3.0° (exp) a b Figure S7: Analysis of replica band intensity in 3.0 • tBG. a Simulated constant energy cut at the Fermi level of 3.0 • tBG overlaid with the mBZ and labels for each mBZ corner. Scale bar is 0.1Å −1 . b Replica intensity of 3.0 • tBG from a, normalised by the intensity of the primary band from the top layer, κ 1 , from experiment and simulation. The contributed intensity comes from the circular areas in a.

EDC analysis of hybridisation gaps
EDCs extracted from cuts along the κ 1 -κ 2 direction were used to determine the size of hybridisation gaps. Examples are shown in Fig. S8 for each of the different twisted graphene stacking arrangements discussed in the main text. Extracted EDCs are fit to a pair of Gaussian functions on a constant background ( Fig. S8d-g). The difference between the Gaussian peak positions is interpreted as the hybridisation gap size, δ.  Figure S9: Simulated and experimental ARPES spectra of 1.5 • tDBG. a-c Experimental (left-hand) and simulated (middle) energy-momentum cuts along the high symmetry directions, as in Fig. 3 of the main text, and the corresponding band dispersions (right-hand). The band dispersions are plotted in the reduced zone scheme of the mBZ. When comparing to the simulated spectra, it is important to recall that the primary bands of each BLG layer are intense across multiple mBZs whilst the intensities of replica bands decrease in higher mBZs and that the broadening of the spectra obscures small gaps and closely spaced bands. The black lines correspond to the peak positions extracted from the experimental data by fitting EDCs and the red lines correspond to the predicted electronic structure. d EDC along the vertical dashed line in the left-hand panel of a. The solid black line is a fit to the data with a pair of Gaussian peaks whose positions correspond to the band positions and separation gives the gap size, ∆ h . e Energy of the simulated (left-hand) and experimental (right-hand) flat-band plotted in the k x − k y plane, with the mBZs overlaid in red. All scale bars are 0.05Å −1 .

Self-consistent analysis of the effect of a back-gate voltage
The effect of a back gate voltage is included in the electronic structure model through a self-consistent analysis that accounts for the change in interlayer potential due to the displacement field and the resultant charge redistribution. The electric displacement field has the following form [4]: where C G is the capacitance to the back gate [17]. This can be used to calculate an initial interlayer potential [18], u, that is the difference between the potential on the top layer, U t , and the bottom layer, U b : where c 0 is the spacing between the layers and ϵ z is the effective out-of-plane dielectric susceptibility. For tMBG, where there are three layers, the energy differences between the two outer layers and the inner layer must be calculated. For this, as in [18], the following parameters are used: c 0,1 = 3.44Å for tBLG, c 0,2 = 3.35Å for BLG, ϵ z,1 = 2.5 for tBLG, ϵ z,2 = 2.6 for BLG. (S16) The initial interlayer potential u i is introduced to the Hamiltonian as follows: (S17) From the wavefunctions calculated using this new Hamiltonian, the new layer density n i in each layer due to the back gate is calculated (n t for the density in the upper monolayer graphene, n m for the density in the middle layer, and n b for the density in the bottom layer of the bilayer graphene), as are the energies of each layer: Here, ϵ hBN = 4 [19], and d = 26 nm is the thickness of the hBN layer for the data in Fig. 5 of the main text. In Eq. (S19), the index i denotes the layer and N is the number of bands being considered for the calculation. Ψ λ are the wavefunctions and f (ϵ l − E F ) is the Fermi distribution. The wavefunctions are used to calculate a new set of interlayer energy differences and the calculations are iterated until they converge and the interlayer potentials are found self-consistently. For tMBG, these give the interlayer potential between the monolayer and the upper layer of the bilayer graphene, u 1 , and between the upper and lower layers of the bilayer graphene, u 2 : Eq. (S19) is then used to calculate the layer densities shown in Fig. 5c of the main text.

Analysis of gated tMBG Dirac cones
To calculate the band parameters as a function of V G from the experimental spectra, as plotted in Fig. 5 of the main text, energy momentum cuts through the κ 1 -κ 2 direction were extracted for each gate voltage (Fig. S10a top panels). Momentum distribution curves (MDCs) were used to extract the band positions of the monolayer and bilayer cones close to the Dirac points. These were fit to Lorentzian functions on a constant background, with the peak centre providing the band position (Fig. S10b). Using standard low energy approximations to the electronic dispersions, the monolayer band positions were fit by E = E D − v|k − k 0 |, and the bilayer band positions by E = E D − 1 2 γ 1 1 + 4v 2 k 2 /γ 2 1 − 1 , Fig. S10a bottom panels, where v is a band velocity. The fitting coefficients provide the Dirac point energy E D . From the tight-binding approximation to the isolated layers, the Dirac point energy can be used to calculate the carrier density, using expressions n MLG = E 2 D πv 2 and n BLG = γ1ED πv 2 for the monolayer and bilayer, respectively, [20]. For simplicity, we have only fit to the valence band, and thus assume the monolayer and bilayer dispersions are symmetric about E D .
The size of the gap at the Dirac point of the bilayer graphene, ∆, is measured in the same way as previously described for the hybridisation gaps, see section S7. An EDC is extracted through the centre of the bilayer cone and fit to a pair of Gaussian functions on a constant background (Fig. S10c). The gap can only be resolved for V G ≥ 7.5 V.